SDE SIR
Simon Frost
using DifferentialEquations
using SimpleDiffEq
using Distributions
using Random
using Plots
using BenchmarkTools
function sir_sde(du,u,p,t)
(S,I,R) = u
(β,γ,δt) = p
N = S+I+R
ifrac = β*I/N*S*δt
rfrac = γ*I*δt
ifrac_noise = sqrt(ifrac)*rand(Normal(0,1))
rfrac_noise = sqrt(rfrac)*rand(Normal(0,1))
@inbounds begin
du[1] = S-(ifrac+ifrac_noise)
du[2] = I+(ifrac+ifrac_noise) - (rfrac + rfrac_noise)
du[3] = R+(rfrac+rfrac_noise)
end
for i in 1:3
if du[i] < 0 du[i]=0 end
end
nothing
end
sir_sde (generic function with 1 method)
δt = 0.01
nsteps = 5000
tf = nsteps*δt
tspan = (0.0,nsteps)
(0.0, 5000)
u0 = [990.0,10.0,0.0]
p = [0.5,0.25,0.01]
3-element Array{Float64,1}:
0.5
0.25
0.01
Random.seed!(1234)
Random.MersenneTwister(UInt32[0x000004d2], Random.DSFMT.DSFMT_state(Int32[-
1393240018, 1073611148, 45497681, 1072875908, 436273599, 1073674613, -20437
16458, 1073445557, -254908435, 1072827086 … -599655111, 1073144102, 36765
5457, 1072985259, -1278750689, 1018350124, -597141475, 249849711, 382, 0]),
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000,
0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x0
0000000000000000000000000000000, 0x00000000000000000000000000000000, 0x0000
0000000000000000000000000000, 0x00000000000000000000000000000000, 0x0000000
0000000000000000000000000, 0x00000000000000000000000000000000, 0x0000000000
0000000000000000000000 … 0x00000000000000000000000000000000, 0x0000000000
0000000000000000000000, 0x00000000000000000000000000000000, 0x0000000000000
0000000000000000000, 0x00000000000000000000000000000000, 0x0000000000000000
0000000000000000, 0x00000000000000000000000000000000, 0x0000000000000000000
0000000000000, 0x00000000000000000000000000000000, 0x0000000000000000000000
0000000000], 1002, 0)
prob_sir_sde = DiscreteProblem((du,u,p,t)->sir_sde(du,u,p,t),u0,tspan,p)
sol_sir_sde = solve(prob_sir_sde,solver=FunctionMap)
retcode: Success
Interpolation: left-endpoint piecewise constant
t: 5001-element Array{Float64,1}:
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
⋮
4992.0
4993.0
4994.0
4995.0
4996.0
4997.0
4998.0
4999.0
5000.0
u: 5001-element Array{Array{Float64,1},1}:
[990.0, 10.0, 0.0]
[989.7575274297931, 10.360050786410875, 0.0]
[989.8182210764063, 10.418767672664574, 0.0]
[989.5704084759752, 10.283557162266963, 0.383023110828716]
[989.3993669198878, 10.47245985961835, 0.36516196956474867]
[989.2332481519006, 10.696048478585872, 0.3076921185844038]
[989.3092616894122, 10.596449490521206, 0.3312775691374906]
[989.2275420377456, 10.35011732186858, 0.6593293894567847]
[989.363627745111, 10.170446479489094, 0.7029145244708828]
[989.3696608412107, 10.080034261645968, 0.7872936462142908]
⋮
[241.03773979754035, 4.480316516535241, 754.8787792103499]
[241.10087333645683, 4.370316992692285, 754.9256451952764]
[241.10586361375567, 4.517953156259023, 754.7730187544108]
[241.06025486079685, 4.399702476740089, 754.9368781868886]
[241.09504471067217, 4.316392938786903, 754.9853978749665]
[241.13547651048168, 4.359867906779514, 754.9014911071644]
[241.15081631338433, 4.370166741338697, 754.8758524697025]
[241.0761466822617, 4.313364127885237, 755.0073247142786]
[240.93476885531518, 4.4336179072290385, 755.0284487618814]
plot(sol_sir_sde)
@benchmark solve(prob_sir_sde,solver=FunctionMap)
BenchmarkTools.Trial:
memory estimate: 670.05 KiB
allocs estimate: 5079
--------------
minimum time: 943.900 μs (0.00% GC)
median time: 1.130 ms (0.00% GC)
mean time: 1.248 ms (4.86% GC)
maximum time: 13.475 ms (88.36% GC)
--------------
samples: 3991
evals/sample: 1